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Abstract 

In  this  paper  we  present  a  method  of  efficiently  implementing  controllers  for  linear  systems  with  large  numbers 
of  sensors  and  actuators.  It  is  well  known  that  singular  value  decomposition  can  be  used  to  diagonalize  any 
real  matrix.  Here,  we  use  orthogonal  transforms  from  the  wavelet  packet  to  “approximate”  SVD  of  the  plant 
matrix.  This  yields  alternate  bases  for  the  input  and  output  vector  which  allow  for  feedback  control  using  local 
information.  This  fact  allows  for  the  efficient  computation  of  a  feedback  control  law  in  the  alternate  bases.  Since 
the  wavelet  packet  transforms  are  also  computationally  efficient,  this  method  provides  a  good  alternative  to  direct 
implementation  of  a  controller  matrix  for  large  systems. 


*This  research  was  supported  in  part  by  a  grant  from  the  National  Science  Foundation’s  Engineering  Research  Centers  Program: 
NSFD  CDR  8803012,  and  by  the  Army  Research  Office  under  the  ODDR&E  MURI97  Program  Grant  No.  DAAG55-97-1-0114  to  the 
Center  for  Dynamics  and  Control  of  Smart  Structures  (through  Harvard  University). 


Figure  1:  Block  diagram  of  control  law  implementation. 


1  Introduction 

This  document  outlines  a  novel  method  of  implementing  controllers  for  dynamic  systems  with  many  (on  the  order 
of  one  thousand)  inputs  and  outputs.  The  idea  presented  here  is  to  perform  a  fast  orthogonal  transform  over  the 
output  space  of  the  system,  changing  the  basis  of  the  output  signal  from  the  Euclidean  basis  (the  standard  basis  for 
II”)  to  the  basis  associated  with  the  transform.  The  hope  is  that  this  new  basis  will  allow  the  implementation  of 
the  desired  control  law  with  fewer  computations  and  less  communication  than  would  be  required  to  implement  the 
equivalent  control  law  in  the  original  basis.  Chou,  Guthart,  and  Flamin  [1]  have  proposed  the  use  of  the  fast  wavelet 
transform  for  this  purpose,  however  these  authors  do  not  suggest  a  systematic  way  of  finding  a  suitable  wavelet  basis. 

Instead  of  limiting  ourselves  to  the  wavelet  basis,  our  approach  uses  an  orthogonal  transform  selected  from  the 
large  collection  of  such  transforms  which  constitute  the  wavelet  packet.  Further,  we  borrow  an  image  processing 
technique  developed  by  Wickerhauser  [4]  to  select  the  “best”  basis  from  the  packet. 

Once  the  basis  is  selected,  the  associated  transform  can  be  used  toward  the  efficient  implementation  of  a  given 
controller.  The  output  vector,  y,  is  transformed  into  the  new  basis  to  yield  the  transformed  output  vector,  y.  The 
transformed  input  vector,  u  is  then  computed  based  on  y.  Finally,  the  actual  control  vector,  v .  is  found  by  performing 
the  inverse  transform  on  u.  This  implementation  is  depicted  in  Figure  1,  where  the  forward  and  inverse  transforms 
are  represented  by  Qi  and  Qo,  respectively. 

The  forward  and  inverse  transforms  from  the  wavelet  packet  are  fast;  an  N  dimensional  vector  can  be  transformed 
in  0(N\ogN )  operations  [3].  In  this  paper,  we  will  show  that  it  is  possible  to  find  a  basis  such  that  the  transformed 
control  vector,  u  =  Ky ,  can  be  computed  in  O(N)  operations.  As  a  result,  the  controller  implementation  depicted  in 
Figure  1  requires  0(N\ogN)  +  0{N )  +  0(N\ogN)  operations.  This  compares  favorably  with  the  0(A7'2)  operations 
required  to  implement  a  controller  directly. 

In  Section  2,  we  introduce  the  plant  model  which  is  used  in  this  paper.  Section  3  describes  the  notion  of  plant 
diagonalization.  It  contains  a  review  of  the  well  known  matrix  factorization  technique  of  singular  value  decomposition 
(SVD)  along  with  a  discussion  of  how  to  approximate  SVD  using  a  wavelet  packet  transform.  The  design  of  a  simple 
controller  which  takes  advantage  of  the  approximately  diagonalized  plant  matrix  is  discussed  in  Section  4.  Finally, 
techniques  presented  in  Sections  3  and  4  are  illustrated  with  an  example  in  Section  5. 

2  Plant  Structure 

For  this  discussion,  we  restrict  ourselves  to  plants  of  the  form 

yc  —  Pcuc  (1) 

where  yc  G  (Cp,  uc  G  <Cm,  and  Pc  G  <Cpxm.  Such  a  system  can  be  used  to  model  any  linear  time  invariant  system 
whose  inputs  are  all  sinusoids  of  a  fixed  frequency,  uj.  Each  element  of  the  input  and  output  vectors  is  a  complex 
number  denoting  the  amplitude  and  phase  of  the  sinusoid.  The  complex  valued  plant  matrix  can  be  thought  of  as  a 
transfer  function  matrix  evaluated  at  the  fixed  frequency  u>. 

The  approximate  diagonalization  algorithm  discussed  in  Section  3  works  only  for  systems  with  real  valued  inputs, 
outputs,  and  plant  matrix  elements.  For  this  reason,  it  is  necessary  to  convert  the  complex  valued  system  of 
Equation  1  to  an  equivalent  real  valued  representation.  This  can  be  done  via  the  following  steps: 

1.  Replace  the  kill  element  of  uc  with  the  column  vector  [real(wc(A;))  imag(wc(fc))]',  for  k  =  1, . . .  ,m. 

2.  Replace  the  A:tli  element  of  yc  with  the  column  vector  [real(yc(k))  imag (yc(k))]\  for  k  =  1, . . .  >p. 

3.  Replace  the  (h,k) th  element  of  Pc  with  the  2x2  matrix 

real  {Pc{h,k))  —imag  (Pc(h,k)) 

imag  (Pc(h,k))  real  {Pc{h,k)) 
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Clearly,  this  transformation  is  invertible.  Inspection  shows  that  complex  multiplication  is  equivalent  to  matrix 
multiplication  under  this  transformation.  Hence,  this  transformation  yields  an  equivalent  real  valued  system  of  the 
form 

V  =  Pu  (2) 

where  y  G  TR2".  u  G  H2m,  and  P  G  IR2px2m. 

Control  design  for  plants  of  this  form  leads  naturally  to  a  “frequency  banded”  approach  to  developing  controllers 
for  more  general  plants.  In  such  an  approach,  the  spectrum  to  be  controlled  is  divided  up  into  frequency  bands  so 
within  each  band  the  plant  can  be  represented  as  a  complex  valued  matrix.  Individual  controllers  are  then  developed 
for  each  band.  These  individual  controllers  can  then  be  added  together  (with  the  proper  filtering)  to  create  a 
controller  which  works  well  over  the  whole  spectrum.  Hence,  the  results  presented  in  this  paper  can  be  applied  to  a 
range  of  plants  which  is  much  broader  than  LTI  systems  with  a  fixed  disturbance  frequency. 

3  Plant  Matrix  Diagonalization 

Recall  that  we  wish  to  find  a  fast  transform  so  that,  in  the  transformed  basis,  the  control  u  can  be  efficiently  computed 
from  the  output  vector  y.  Singular  value  decomposition  (SVD)  can  be  used  to  generate  a  good  basis  for  our  purposes. 
SVD  can  be  used  to  diagonalize  any  real  matrix  [2].  Consider  the  real  valued  p  x  m  matrix  P.  SVD  states  that  P 
can  be  factored  such  that 

P  =  Q^Ql,  (3) 

where 

1.  The  columns  of  Q\  G  IRpxp  are  the  eigenvectors  of  PPT . 

2.  The  columns  of  Qo  G  IRmxm  are  the  eigenvectors  of  PTP. 

3.  The  p  x  to  matrix  S  is  diagonal,  and  values  on  the  diagonal  are  the  square  roots  of  the  eigenvalues  of  PPT 
and  PTP. 

The  matrices  Q\  and  Q 2  are  orthogonal,  i.e.  QiQf  =  Ip  and  Q2Q2  =  Im,  where  In  is  the  n  x  n  identity  matrix. 
Hence,  multiplying  P  by  Qf  on  the  left  and  Qo  on  the  right  yields 

Q?>Q2  =  E.  (4) 

In  the  present  context,  P  is  a  plant  matrix  of  an  m-input,  p-output  system.  In  other  words, 

V  =  Pu,  (5) 

where  u  G  Hm  is  the  input  vector  and  y  G  IRP  is  the  output  vector. 

Now  let  y  =  Qjy  and  u  =  Qju.  Substitution  into  the  plant  equation  (Equation  5)  yields 

y  =  Q  I  PQ211 

=  m.  (6) 

The  ith  element  of  the  transformed  input  vector  u  affects  only  the  ith  element  of  the  transformed  output  vector  y 
for  i  =  1,2,..., min(p, to).  The  entire  system  can  be  controlled  using  local  feedback,  i.e.  each  input  is  determined 
based  solely  on  the  value  of  its  corresponding  output.  In  the  transformed  basis,  the  control  u  can  be  computed  from 
y  in  0(min(p, to))  operations. 

Unfortunately,  the  SVD  generated  transforms  are  computationally  intensive;  it  takes  0{N'2)  operations  to  trans¬ 
form  an  N  dimensional  vector.  This  obstacle  can  be  overcome  by  using  a  fast  wavelet  packet  transform  (FWPT) 
which  approximates  the  Karhunen-Loeve  transform.  As  will  be  shown  later  in  this  section,  the  Karhunen-Loeve 
transform  can  be  used  to  perform  SVD.  An  algorithm  to  find  such  a  transform  has  been  presented  by  Wickerhauser 
[4].  The  FWPT  is  computationally  efficient;  it  takes  0(N\ogN)  operations  to  transform  an  N  dimensional  vector. 

We  know  that  the  Karhunen-Loeve  representation  of  a  given  ensemble  of  vectors  has  a  lower  entropy  than  any 
other  orthonormal  representation.  Wickerhauser ’s  algorithm  recursively  searches  all  of  the  orthogonal  bases  in  the 
wavelet  packet  to  find  the  basis  in  which  the  vector  ensemble  has  the  lowest  entropy  representation.  This  basis  is 
then  “closer”  to  Karhunen-Loeve  than  any  other  basis  in  the  packet.  Hence,  the  wavelet  packet  transform  to  this 
basis  will  be  used  as  our  “approximate”  Karhunen-Loeve  transform. 

If  we  consider  an  ensemble  of  vectors  composed  of  the  columns  of  P,  the  Karhunen-Loeve  transform  for  this 
ensemble  is  identical  to  the  SVD  factor  Qf.  Similarly,  if  we  let  the  ensemble  be  the  rows  of  P,  the  resulting 
Karhunen-Loeve  transform  is  identical  to  the  SVD  factor  Q.[  .  Thus,  we  can  use  the  Karhunen-Loeve  transform  to 
diagonalize  the  plant  matrix  P.  Or  we  can  use  an  “approximate”  Karhunen-Loeve  transform  to  “approximately” 
diagonalize  P. 
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Figure  2:  Block  diagram  of  full  feedback  system. 


4  Controller  Design 


Before  we  discuss  the  development  of  the  feedback  control  law  for  this  method,  we  review  the  centralized  case  so  we 
have  something  to  which  we  can  compare  the  results.  Consider  the  feedback  system  depicted  in  Figure  2.  Here,  d 
is  the  disturbance  input,  u  is  the  control  input,  and  y  is  the  output.  The  object  is  to  find  a  control  matrix  K  such 
that  the  effect  of  the  disturbance  on  the  output  it  minimized.  Consider  the  case  where  P  is  square  (i.e.  p  =  m)  and 
invertible.  Then  letting  K  =  7P-1  yields 

V  = 


This  control  law  reduces  the  transmission  from  d  to  y  by  a  factor  of  1  +  7. 

Now  consider  the  system  shown  in  Figure  1.  If  we  let  E  be  the  approximately  diagonalized  plant  matrix,  then 

simple  matrix  manipulation  shows  that  =  7-P-1-  As  a  result,  letting  K  =  7E-1  accomplishes  a 

disturbance  transmission  reduction  by  a  factor  of  1  +  7. 

Unfortunately,  E  is  only  “approximately”  diagonal,  so  7E-1  will  not  be  a  diagonal  matrix.  This  fact  ruins 
the  locality  of  the  feedback  law  in  the  transformed  space.  This  locality  is  the  feature  which  allows  the  efficient 
computation  of  u.  Hence,  we  need  to  find  a  way  to  implement  something  close  to  7E-1  while  preserving  as  much 
locality  as  possible. 

The  first  idea  that  comes  to  mind  is  to  simply  ignore  the  off-diagonal  elements  of  E.  In  other  words,  define  E 
to  be  equal  E  with  all  off-diagonal  elements  replaced  by  zeros.  Now  the  feedback  law  7E-1  is  diagonal  and  can  be 
implemented  locally.  In  sequel,  we  will  denote  this  feedback  law  as  the  “pure  local  controller” . 

—  T  — 

The  pure  local  feedback  law  will  work  well  if  the  transforms  Q 1  and  Q 2  are  good  approximations  of  Q{  and 
Q-2,  but  this  is  not  always  the  case.  If  the  approximate  transforms  are  not  sufficiently  close  to  the  desired  Karhunen- 
Loeve  transforms,  then  E  will  have  off-diagonal  elements  which  cannot  simply  be  ignored.  In  this  case,  we  will  be 
forced  to  violate  the  locality  of  the  feedback  by  allowing  a  small  number  of  the  transformed  inputs  have  access  to 
the  information  from  a  small  number  of  the  transformed  outputs. 

The  first  step  in  doing  this  is  to  set  all  of  the  off-diagonal  elements  of  E  whose  absolute  values  are  below  a  certain 
threshold  to  zero.  Let  n  be  the  number  of  non-zero  off-diagonal  elements  which  remain.  Next,  the  transformed 
inputs  and  outputs  are  reordered  so  that  the  resulting  matrix  takes  the  block  diagonal  form 


P(d  -  7 P^y) 
P  d  -  yy 


Mi  0 


E  = 


0 


(7) 


where  M,;,  i  =  1, . . . ,  k,  k  <  n  are  square  and  S  is  (m  —  n)  x  (m  —  n )  and  diagonal.  The  inverse  of  E  is  then 


E-1 


M, 


-1 


0 


0 


(8) 


As  a  result,  to  implement  the  control  law  K  =  7E  1 ,  a  few  of  the  inputs  require  information  from  a  few  of  outputs 
but  the  rest  of  the  inputs  can  be  calculated  locally.  We  will  denote  this  feedback  law  as  the  “local-Hr  controller” . 
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Figure  3:  Diagram  of  4  x  4  rectangular  array  of  CSM  Tiles. 

The  hope  is  that  feedback  of  the  pure  local  or  local+n  (with  n  small)  form  will  perform  comparably  with  the  full 
feedback  law.  For  this  comparison,  we  look  to  the  example  in  the  next  section. 

5  Example 

The  Composite  Smart  Material  (CSM)  Tile  [5]  provides  the  motivation  for  the  following  example.  The  CSM  tile 
is  a  self  contained  unit  containing  underwater  acoustic  sensors  and  actuators  along  with  a  power  supply  and  signal 
processing  hardware.  In  practice,  thousands  of  these  “smart  tiles”  will  be  mounted  on  the  outside  of  the  hull  of  a 
submarine,  covering  the  entire  hull.  This  massive  sensor/actuator  array  will  then  be  used  to  actively  reduce  acoustic 
radiation,  cancel  enemy  sonar  pulses,  and  perform  acoustic  sensing. 

Here  we  consider  the  example  of  a  rectangular  4x4  array  of  the  CSM  Tiles.  We  assume  that  the  disturbance  signal 
is  a  sinusoid  of  a  fixed  frequency,  d  =sin(27r /),  so  that  we  can  model  that  interaction  between  the  tiles  as  a  complex 
valued  plant  matrix.  The  medium  between  the  tiles  is  homogeneous  so  that  the  wavelength  of  the  disturbance  signal, 
A,  is  fixed.  The  magnitude  of  the  coupling  between  any  two  tiles  is  where  d  is  the  distance  between  them.  The 
phase  shift  between  any  two  tiles  is  27T  (d  mod  A).  These  matrices  are  shown  in  the  top  half  of  Figure  4  where  the 
real  valued  elements  of  the  matrix  are  replaced  by  pixels  of  varying  intensity. 

The  first  step  is  to  transform  the  complex  valued  plant  matrix  to  its  real  matrix  equivalent  as  described  in 
Section  2.  This  yields  a  32  x  32  real  valued  matrix.  The  absolute  values  of  the  elements  of  this  matrix  are  shown 
in  the  bottom  half  of  Figure  4.  Now  we  can  write  a  model  of  the  system  in  the  form  of  Equation  2.  The  input 
vector  u  and  the  output  vector  y  are  both  in  IR32.  The  real  and  imaginary  parts  of  the  sensor  output  of  the  A:tli  tile 
become  y{2k  —  1)  and  y(2k),  respectively.  Likewise,  the  real  and  imaginary  parts  of  the  input  to  the  actuator  on 
the  fcth  tile  become  u(2k  —  1)  and  u(2k),  respectively.  This  means  something  important  to  the  local+n  controller 
design:  since  each  sensor  effectively  measures  two  outputs  (real  and  imaginary  part)  and  each  actuator  effectively 
takes  two  control  signals,  adding  one  connection  between  two  tiles  effectively  brings  the  data  from  2  inputs  and  2 
outputs  together.  This  means  that  the  local+n  controller  will  include  the  2x2  blocks  with  the  largest  norm  in  the 
approximately  diagonalized  matrix. 

The  next  step  is  to  find  the  approximate  forward  and  inverse  transforms  to  diagonalize  the  plant  matrix  P.  This 

was  done  and  the  transforms  were  used  to  compute  the  approximately  diagonalized  plant  matrix,  S  =  Qi  PQ o. 
The  absolute  values  of  the  elements  of  this  matrix  are  shown  in  Figure  5. 

Controllers  for  this  system  were  developed  as  outlined  in  Section  4.  The  parameter  7  was  set  to  9  so  that  the 
full  feedback  control  case  attenuates  the  disturbance  by  a  factor  of  10.  The  disturbance  transmission  matrices  for 
the  local+5  controller  are  shown  in  Figure  6.  For  this  controller,  Figure  5  was  used  to  determine  which  off-diagonal 
blocks  to  include.  As  the  figures  show,  the  performance  of  the  local+5  controller  is  very  close  to  that  of  the  full 
feedback  control. 

Intuition  tells  us  that  as  n  gets  larger,  the  performance  of  the  local+n  controller  should  get  better.  To  check 
this,  we  introduce  a  measure  of  closeness  of  the  performance  of  the  local+n  controller  to  the  performance  of  the  full 
feedback  controller.  Let  Pfuii  and  Pn  be  the  closed  loop  disturbance  transmission  matrix  under  full  feedback  control 
and  local+n  control,  respectively.  We  now  define  the  error  of  the  local+n  controller,  En,  as 

En  =  \\Pfuii  ~  Pn\\2  (9) 

where  ||-||2  represents  the  matrix  2-norm,  or  largest  singular  value.  A  plot  of  En  vs.  n  is  shown  in  Figure  7.  The 
plot  matches  our  intuition,  but  it  also  shows  something  that  is  not  obvious.  In  a  16-input,  16  output  system,  there 
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Figure  4:  Plant  matrices  for  4  x  4  rectangular  CSM  Tile  array. 
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Figure  5:  Approximately  diagonalized  plant  matrix,  E  =  Qi  PQ2 . 
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Figure  6:  Disturbance  transmission  matrices  for  local+5  controller. 
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Figure  7:  Error  of  local+n  controller,  En  =  \\Pfuii  —  Pn||0 


are  (16  x  16  —  16)/2  =  120  possible  off-diagonal  connections  to  be  made.  In  other  words,  the  local+120  controller  is 
equivalent  to  the  full  feedback  controller.  But  as  Figure  7  shows,  the  error  of  the  local+n  controller  is  zero  for  all 
n  >  56.  This  means  that  for  this  example,  it  is  possible  to  exactly  reproduce  the  performance  of  the  full  feedback 
controller  using  less  than  half  of  the  possible  connections. 

6  Conclusion 

Here  we  have  presented  a  computationally  efficient  method  of  implementing  controllers  for  large  scale  systems.  In  this 
method,  the  output  vector  is  transformed  into  a  basis  which  allows  full  feedback  control  using  only  local  information. 
This  can  be  done  because  the  basis  transform  is  chosen  to  approximately  diagonalize  the  plant  matrix.  The  resulting 
locality  provides  advantages  in  both  design  and  implementation.  These  results  should  prove  useful  for  the  growing 
number  of  applications  which  require  large  numbers  of  sensors  and  actuators. 
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